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Abstract 



In this paper, we propose a numerical regularized moment method to solve the 
Boltzmann equation with ES-BGK collision term to simulate polyatomic gas flows. 
This method is an extension to the polyatomic case of the method proposed in [5], 
which is abbreviated as the NRa:a; method in [S, • Based on the form of the Maxwellian, 
the Laguerre polynomials of the internal energy parameter are used in the series ex- 

■ pansion of the distribution function. We develop for polyatomic gases all the essential 
techniques needed in the NRxx method, including the efficient projection algorithm 

■ used in the numerical flux calculation, the regularization based on the Maxwellian it- 
"j^ , eration and the order of magnitude method, and the linearization of the regularization 

term for convenient numerical implementation. Meanwhile, the particular integrator 
in time for the ES-BGK collision term is put forward. The shock tube simulations 
with Knudsen numbers from 0.05 up to 5 are presented to demonstrate the validity 
, , of our method. Moreover, the nitrogen shock structure problem is included in our 

^ • numerical experiments for Mach numbers from 1.53 to 6.1. 

^3 ■ Keywords: polyatomic ES-BGK model; moment method; NRxa; method 

§ ■ 1 Introduction 
(N" 

The kinetic theory has long been playing an important role in the rarefied gas dy- 
namics. As a mesocopic theory standing between the fluid dynamics and the molecular 
^ . dynamics, the kinetic theory is built on the basis of the Boltzmann equation, which uses 

I a distribution function to give a statistical description of the distribution of microscopic 

particles' velocities. In 1940s, Grad [TT] proposed the idea using the Hermite expansion to 
approximate the distribution function, and a 13- moment theory was given in detail in jllj . 
Recently, based on the idea of Grad, systems with large numbers of moments together with 
their numerical schemes are considered in O [9l [8] , where some regularizations inspired 
by [21\ 119] are also taken into account. In [H], the numerical regularized moment method 
is abbreviated as the NRxx method. However, all these works concentrate only on the 
monatomic gases, and in this paper, we will develop the NRxx method for the polyatomic 
case. 

The study to apply the moment method to polyatomic gases can be traced back to 
McCormack [ISj, where a 17-moment model was proposed. As far as we know, the most 
recent polyatomic extension of Grad's 13-moment equations is the work of Mallinger [1] 
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whose system contains only 14 moments. In both models, a great amount of work is 
devoted to the deduction of the collision terms. In order to generalize the moment theory 
to large number of moments, we prefer a BGK-like simplified collision operator. As in the 
monatomic case, the simplest BGK model fails to give correct heat conduction, and for 
polyatomic gases, it also gives incorrect relaxation collision number, resulting in qualitative 
errors in temperatures compared with the Boltzmann equation [3]. Possible alternatives 
include the Rykov model [18] and the ES-BGK model [3], which incorporate physical 
Prandtl number and relaxation collision number into the collision term. In this work, our 
investigation is restricted to the ES-BGK model. 

For polyatomic gases, besides the velocities of microscopic particles, an additional 
nonnegative ordinate representing the energy of internal degrees of freedom appears in 
the distribution function. Thus, in order to expand the distribution function into series, 
the basis functions are chosen as a combination of Hermite polynomials and Laguerre 
polynomials with proper translation and scaling based on the macroscopic velocity and 
translational and rotational temperatures of the gas. By considering the coefficients of 
the basis functions as moments, a system with infinite number of moment equations is 
derived. A moment closure is then followed to truncate the system with infinite equations 
and get a system with only finite equations. The framework for the moment closure is the 
same as [9]: 

1. the Maxwellian iteration is applied to determine the order of magnitude for each 
moment; 

2. by dropping higher order terms, the truncated moments are expressed by moments 
with lower orders; 

3. for easier numerical implementation, the expression is linearized around a Maxwellian. 

However, the details of the Maxwellian iteration are significantly different. In the poly- 
atomic case, the iteration is much more complicated than the monatomic case because of 
the existence of both translational and rotational temperatures in the basis functions, and 
the process should be conducted carefully. Moreover, for the ES-BGK model, analysis 
on the moments of the Gauss distribution also increases the complexity. Fortunately, the 
final result remains a similar form as simple as in [9]. 

As to the numerical scheme, the general framework in [5] is applicable. A split scheme is 
applied to divide the transportation part and the collision part, and the transportation part 
is processed by a finite volume method. Recalling that a special "projection" introduced in 
[3 [8] is required in the calculation of numerical fluxes, we further develop this technique to 
the polyatomic case in this paper. Meanwhile, the polyatomic ES-BGK collision term can 
no longer be solved analytically as the BGK operator [7], and the Crank-Nicolson scheme 
is applied to ensure the unconditional numerical stability. Our numerical experiments 
show that our scheme correctly converges to the solution of the Boltzmann equation as 
the number of moments increases. The distinction between BGK and ES-BGK models, 
together with the relation between monatomic and polyatomic cases, is illustrated by the 
numerical results of shock tube problems. Also, we apply the NRxx method to the nitrogen 
shock structure problem, and the results are comparable to the experimental data. 

The rest of this paper is arranged as follows: in Section [21 a brief review of the poly- 
atomic ES-BGK model is given. In Section [31 the polyatomic NRxx method is introduced 
comprehensively, and in Section [H a number of numerical experiments are carried out to 
validate our algorithm. As a summation, some concluding remarks are given in Section [5l 
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Finally, some involved calculations are collected in the appendix for better readability to 
the body matter. 



2 The ES-BGK Boltzmann equation for polyatomic gases 

The ES-BGK model for polyatomic gases, which gives correct Navier-Stokes heat con- 
duction compared with the BGK model, has been deduced in [6]. The polyatomic 
ES-BGK Boltzmann equation reads 

^+^-V,/ = Pr.^(G-/), (2.1) 

where / denotes the molecule distribution, which is a positive function / = f{t,x,$,,I) 
with a;, ^ E and t, / G M"*". The parameters t, x and ^ stand for the time, spatial position 
and microscopic molecule velocity respectively, and / is an internal energy parameter. In 
the right hand side of ()2.ip . Pr is the Prandtl number, p is the pressure, and /i denotes 
the viscosity coefficient. G is a generalized Gaussian defined as 

Here 5 is the total number of molecular internal degrees of freedom, and R is the gas 
constant. The density p and the macroscopic velocity u are related to the distribution 
function / through 

p = [ /d^dl, u=- [ Ud$dl, (2.3) 



and Tj-ei is a relaxation temperature 

T,el = Z-^Teq + (1 - Z-^)Ti,,t, (2.4) 

where Z is the relaxation collision number. For polyatomic gases, three temperatures are 
used frequently, including the translational temperature Ttr, the internal temperature Tint, 
and the equilibrium temperature Tgq. They are defined by 

Ttr = ^[ \$-u\^fd$dI, (2.5) 
^-t = A/ ^'^'/d^dl, (2.6) 

Teq = (3rto + 5Tint)/(3 + 5). (2.7) 

And the pressure p is obtained from the ideal gas law: 

p = pRTeq. (2.8) 
Now it only remains to define and T: 



where 



/ 



r = (1 - Z-^)[{1 - iy)RTtrId + uQ/p] + Z-^i^Teqld, (2.9) 



e=/ {i-u)(^{^-u)fd^dl, u = \ ^J_. , (2.10) 



and Id stands for the identity matrix. 
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3 The NRcca? method for polyatomic ES-BGK equation 



In this section, we are going to extend the NRxx method proposed in ^ i9j to the 
polyatomic case, which includes the following steps: 

1 . The distribution function is expanded into a series with specially selected basis func- 
tions. 

2. A system with infinite number of moment equations is deduced. 

3. The moment system is truncated at a certain place and made closed by regulariza- 
tion. 

4. The regularization term is linearized in order to simplify the numerical implementa- 
tion. 

5. The numerical method is carried out following [5]. 
The details are introduced in the following five subsections. 

3.1 Spectral representation of the velocity space 

In the NRxx method for the monatomic gases, the Hermite polynomials have been 
employed to construct the basis functions of the velocity space, since Hermite polynomials 
are orthogonal over the region (— oo, +oo). For the polyatomic distribution function, since 
/ € M'*', we use the Laguerre polynomials, which are orthogonal over the region [0,+oo), 
as the basis functions in the ordinate /. Thus the basis function has the following form: 

V.,.,T.,T,„,(^, J) = ] (7f^)"'(i?Ti.t)-(^/^+'=)4"^(J)exp(-J) . 

(V2^j {RTt,)-— n Heo^A^d) exp ( 

d=l ^ ^ 

where a = (ai, a2, as) is a multi-index, and 

T-m„J jfc 

Lt\j) = -J^fj,ie-'j''-n, (3.3) 
He^ix) = {-ir exp (:^) exp 1^-^^ . (3.4) 

Some properties of the Laguerre polynomials L^j^^ and the Hermite polynomials Hen can 
be found in AppendixEl With equation ()3.ip . the distribution function /(^, /) is expanded 
as 



(3.5) 

Let us consider the general case when T^^ , Tint and u have no relation with the distribution 
function /. Using the orthogonality of the Laguerre and Hermite polynomials, we can 
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deduce that 

/d^d/ = /o,o, (3.6a) 



/ 



f 



e,/d^d/ = /o,on, + /e,,o, i = 1,2,3, (3.6b) 



/ d^ dl = IdRTi^tfo.o - hi, 



(3.6c) 



/ 



h^l^fd^dl = ^Ofiluf + (^RTtrf0,0+Ujfe„0 + ^e^o) • (3.6d) 



If u is the macroscopic velocity and Ttr, Tint are the translational and internal temperatures 
for the distribution /, using ()2.3p . ()2.5p . ()2.6p and ()3.6p . we conclude 

3 

/0,0 = /5, /ej ,0 = /0,1 = ^ /2ed,0 = 0, j = 1, 2, 3. (3.7) 

d=l 

If (|3.7p is satisfied, then ()3.5p is called as a normal representation of f. If ()3.5p is not a 
normal representation, then the density, momentum and translational and internal energies 
can be easily calculated through (|3.6p . For a normal representation, we have 

Qij = (1 + 5ij)fe^+ej,o + 5ijpRTtr, ij = 1, 2, 3. (3.8) 
where is defined in ()2.10p . 

3.2 The moment equations for the ES-BGK model 

In this section, we are going to derive equations for the moment set {fa,k}- The general 
strategy is to substitute (j3.5p into (j2.ip , and then match the coefficients of the same basis 
functions. For the left hand side of ()3.5p . the process is similar as that in [9j, and the 
detailed derivation can be found in Appendix [Bl Suppose G has the following expansion: 

Then the analytical expressions of the moment equations are obtained as 

3 



(3.9) 



dfa,k ,^dUd. , ld{RTtr)^. . , ,,a(iirint) „ 

+ 2^ -gfJo,-ea,k + 2 ol 2^ Ja-2ea,k - [m + k) — ^,fc-l 

d=l d=l 

E\ f ^h-ej,k dfa^k , , , sdfa+ej,k\ 



i=i 

+ ^ 'Q^. {RTtrfa-ea-ej,k + Ujfa-ea,k + {aj + I) fa-e^+ej ,k) 

+ 2~^Q~~^ ^ {RTtrfa-2ea-ej,k + Ujfa-2ea,k + (Oj + l)/a-2ed+ej ,fc) 

- ("I + ^)~^^^ {RTtrfa^e,,k^l + Ujfa,k-1 + ("j + fa+e, ,k-l) 

= Pr ■ —{Ga.k ~ fa.k), 
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where /^^; is taken as zero when / or any of the components of /? is negative, and m is 
defined in 

Now we focus on the relation between Ga,k and fa,k- The expression of G ()2.2p — 
(j2.10p and the equahties under normal representation (j3.7p and (|3.8p show that Ga,k are 
functions of p, Ttr, Tint and /e^+ej,o with i,j = 1, 2, 3. Thus the system ()3.10p is closed for 
a G and k = 0,1, which means only the expressions of Gafl and Ga,i are needed. The 
following results are trivial: 



Go,o= / Gd^dl = p, (3.11a) 

Ge„o= [ CjGd$dI-Go,oUj = 0, (3.11b) 
Go,i = ^-RTintGofi - [ /2/^Gd^ dl 

= ^-p\RT,r.t - Z-^RT,^ - (1 - Z-i)i?rint] (3-llc) 



2Z 



p{RTmt — RTeq), 



where ()3.11cp comes from the physical meaning of the relaxation collision number. More- 
over, since G(^, /) has the form 

G{CI) = Gi{$)G2iI), (3.12) 

where 



(3.13) 

and the basis functions a, k, Ttr, Tint have the same structure (see (jB.ip ) 

i^a,k,Ttt,Tinti'"^ J) = ^l,a,TtAv)^p2,k,TintiJ)^ (3-14) 

we can conclude 

Ga,i = CGa,o, (3.15) 
where C is independent of a. From (|3.11ap and ()3.11cp . we find 

C= — (flTi„t-i?Teq). (3.16) 

Until now, what remains is to work out the expressions of Ga,o for \a\ ^ 2. This needs 
some involved calculation with details in Appendix O The final result is in a recursive 
form as 



1 ^ 

Ga,o = - Y^[{1 - Pr-i)(ei,/p - RTtAj) + Z-\RT,q - i?Ttr)<5,j]G„_e,-e,,o, (3.17) 



where i G {1,2,3} such that Qj > 0, and Ga-e^-e ,o is taken as zero when Uj — 5ij — 1 < 0. 
With ()3.11bp . one can easily observe that Gafl = when |q| is odd. 
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As the end of this subsection, we give the equations for the velocity u and the trans- 
lational and internal temperatures Ttr, Tint- The equation for can be obtained by 
substituting a with in (|3.10p . and the result is 



P^ + L.[P^^9^^+^)=^^ ^ = 1'2,3. (3.18) 

The equation for Ttr can be obtained by substituting a with 2ei, 2e2, 2e3, and then 
summing up all the three equations. The result is 

where 

3 

=2/3e,,0 + J^/e,+2e„0, J = 1, 2, 3. (3.20) 

d=l 

Similarly, if we set a = and k = 1, then we have 

1=1 •' 7=1 ■' 



3.3 Truncation and closure with regularization 

Since the moment system (j3.10p contains an infinite number of equations and cannot be 
used for computation, we need to choose a finite set from them as the governing equations 
of our method. However, due to the existence of the last term in the second line of (|3.10p . 
the resulting moment system will be unclosed, which leads to the "closure problem" for 
the moment method. 

We first consider the truncation of the spectral expansion. In general, we can choose 
two non-negative integers Mq ^ 2 and Mi ^ 0, and use the moment set {/a,o}|a|^Mo 
{fa,i}\a\^Mi the finite subset. Such choice well retains the Galilean invariance since 
{fafl} and only couple with each other in the collision term. We will postpone the 

discussion of the relation between Mq and Mi. Below we use X to denote the index set 
such that the set {fa,k}{a,k)el contains all the moments appearing in the final moment 
equations. 

Now we are going to make the system closed. The simplest way is to set fa+e ,k = 
in (j3.10p if (a + ej,k) ^ I, which leads to the Grad-type moment equations for poly- 
atomic gases. Mallinger's work }14j has generalize the Grad 13-moment equations to the 
polyatomic case. Here, we follow [n\ \T9\ [9] and use the Maxwellian iteration together 
with the order of magnitude method to give a more reasonable closure. The procedure of 
Maxwellian iteration is constructed by rearranging ()3.10p and adding superscripts repre- 
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senting the iteration steps: 



a Ai^) 3 n 
'^Jg^k OUd An) 

dt ^ dt ■'"-^d.fc 



d=l 



+ 



2 dt 

3 



d=l 



tr 



9xj 



9Xj 



I5(i?7]^^^ A / („+!) („) („) inf(") 

~'~ 9 /^T- . I -fi-^tr Ja-2ea-ej,k "■jJa-2ea,k "i" V«i "r ^;7a_2ed+ej 



9Xj 



where e = /x/(Pr -p) is considered as a small parameter and {a,k) satisfies 
(a, A;) G 5 := {N^ x {0, 1} : |a| ^ 2 if /c = 0, |a| ^ 1 if = 1}. 
In the first line of ([322]), a'-^I and 5(e) are defined as 



A 



H 

a,k 



^Pi ■ Z p{RTcq — -RT^j, if a = e, + Cj, k = 0, 

other cases. 



'^a,k^ 



B{e) 



Pr • e if \a\ =2, k = 0, 
£, other cases. 



(3.22) 
(3.23) 

(3.24) 
(3.25) 



The reason why the iteration scheme for \a\ = 2 is special is that Ga,k are fmictions of 
fa, |a| = 2. Note that /o,o(/o)i u, /e^,^o(= 0) and /o,i(= 0) remain invariant during the 
iteration, and according to ()3.19p and ()3.2ip . Ttr and Tint evolves as follows: 



^iLr'^ = Tec, - eZ 



dii^ A dT^ 2 ^ 

i=i ^ ^ .7=1 



Ox 



1 .7 = 1 



where Teq also remains invariant during the iteration. The iteration starts with 



T,?) = = r,q, Z^;;^ = p, = for (a, fc) / (0, 0). 

In order to simplify the notation, we define the following vectors: 



f{0) 



.(0) 



m,k 



\J a,kJ\c(\=rm ^ [. 



(n) 

mi ,r?i2] ,k 



An) 



mi < 7712. 



(3.26) 
(3.27) 

(3.28) 
(3.29) 
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and rewrite ()3.22p and ()3.26p as 

Ja,k — ^a,k ^'-a,k ' [|a|-3,|Q!|+l],fe ^-^a.fc ' [|a|-l,|Q:|+l],fc-l' 1,0. OU; 

rp(n+l)_rp _rr(T'-"-'l ^ Ci-i-W 

^tr — ^eq ^'^K^tr ' "'^ [2,3] ,0'' ' {O.Ol) 

Tt^'^=T.q-eC{Tt^,Ftl), (3.32) 

where £(•, •) is a linear operator, and c}~^\ is a vector of linear operators. Each of ;C^"^'s 
components has the following form: 

E E T.c, J ' l.,^ . -{o,l)^ (3.33) 

where Cj^r,s = Cj^r,sip,u,d/dt,Vx), which can be considered as a "constant" during the 
iteration. L^^j^ is a vector, whose components can be expressed as 

E E cJ-^^^[RT^^r%\ ^e{o,ir, (3.34) 

j = l Si+S2 = l, 53+54^1 ^ 

where Cj^s are also constants. Now we are ready to carry out the Maxwellian iteration. 

The first step of iteration In the first step, the formulae for the translational and 
internal temperatures can be written as 

r,(;) = Teq + e5(°) , T^, = Teq + ^C/^") , (3.35) 
where 5^°) ~ 0(1) and C/(°) ~ 0(1). It is easy to find 

Gi°i = /i% VaGN^, A: = 0,1. (3.36) 

Thus 

^{0) ^ ( -lePv-Z'^pRS^-'^l ifa = 2ei,k = 0, 
"''^ \ 0, other cases. 

Now, it can be easily deduced from (|3.30p that fj^l. is nonzero if and only if G [|ce| — 
3, \a\ + 1], A; = or G [\a\ — 1, \a\ + 1], k = 1. Precisely, we have 



0(e), (a, k) £ S and \a\ ^ 3, A; = 0, 
0(e), {a,k) G 5 and \a\ = 1, k - 
0, other cases for (a, k) G 5, 



fal ~ { 0(e), (a, k)eS and |q| = 1, = 1, (3.38) 



and 

[2,3],0 ~ [2,3],0 "^^-"[2,3],0' 1,1 " ^ 1,1 1,1, 1,0. OJJ 

where -^^[2 3] k -^1*^1 order of magnitude 0(1). The meaning of the subscripts 



Now let us consider G^^\. Go,o> Go,i and Ga^k with odd \a\ keep invariant during the 



[2 _ 

of H is the same as F. 

Now let us consider 
iteration. For Gafi with \a\ ^ 2, p.l7p gives 

3 

^(n) _sr^ J (rpin) rp -p{n)\ f^{n) (1 An\ 

^afi - 2^ yj^tr " -^eq, ^ 2,0 ) a-a-Sj ,0^ 

3=1 
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where La,j{-, •) is a linear function independent of e. Using ()3.38p . ()3.35p and Go,o = P 
0(1), we can easily get 

Gafi ~ 0(e'"'/^), |a| is even. 

Using ()3.16p . we have 



\a\ IS even. 



(3.41) 
(3.42) 



The general progress In general case, we have 



An+l) 
J a.k 



~ < 



' 0(en"l/3l )^ (q,^ k) eS and \a\ ^ 3(n + 1), A; = 0, 

o(en"l/3l )^ la, k) gS and \a\ >3{n + l),k = 0, 

0(er(l"l+2)/3l )^ k)eS and \a\ ^ 3n + 1, A: = 1, 

o(er(l"l+2)/3l )^ (q,^ A;) G 5 and \a\ >3n + l, k = l. 



(3.43) 



and 



-,(n+l) _ „{n) n+1 rrW 

[2,3(n+l)],0 ^ [2,3(n+l)],0 ^ [2,3(n+l)],0' 

p(n+l) _ p(n) n+1 tr{n) 

[l,3n+l],l [l,3n+l],l ^ [l,3n+l],l' 



tr 

(n+1) 



mt 



'tr 
' int 



(3.44) 
(3.45) 
(3.46) 
(3.47) 



where h''^^ ^ ~ 0(1), ~ 0(1) and ~ 0(1). (IMjl — (13:^71) can be validated 
by induction. Through the derivation in the last paragraph, we have known that ()3.43p — 
(IHITT)) hold for n = 0. Now we suppose (fTIHD — (fXITD hold forn - 1. Then (fXIOD gives 



G 



3 = E[- 



i=i 



-teq, 2 



~r fc -'^aj I '-"20 



(n) 



(3.48) 



Now using Gq"q = G, 



.(n-l) 
0,0 



p and ()3.40p with n replaced by n — 1, simple induction gives 



^(n) _ ^(n-1) 
^a,0 "-^0,0 



0{e 



n+\a\/2~U 



\a\ ^ 2 and IqI is even. 



(3.49) 



For G^^\, we have 



(-,(n-l) , ^p.nrr{n-l) 

^ 2Z 



^(n-l)^(n-l) 



^(n-l)(c(n) _ ^(n-D) ^ ^^e^U^-'^ G^:l 

0{e) ■ 0(e"+l"l/2-i) + . 0(1) • 0(el"l/2) 
0(e"+l°l/2), |q| ^ 2 and \a\ is even. 



Similarly, by (|3.3ip . we have 



(n+1) 



tr 



p (rp{n-l) Tp(n-l)\ , n r ( q(n-l) rrin-l) 
I ^tr '-^ [2,3],0 i ^ '-y^ '-"[2,31,0 



rp{n) _ n+1/- ( Q(n-l) rj-(n-l) 
-^tr £ I '-"[2,3],0 



(3.50) 



(3.51) 
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Defining = -£ (^("-i), //[^ 3]^^) gives (IMjl . The validation of i^^:^ is almost the 



same. I^M), and (l33T]l imply that A^"^ 

Now it remains to consider ^ 
(I33T]) show 

(n+l) ^(n) /ore-^+l^ 
(n+1) j-{n) 

Using the assumptions of the induction, one has 



never greater than 0{e'' 



a,k 

(n+1)^ Similar as p.SOp . the equations (|3.33p . ()3.34p and 



C 
L 



-(n+1) 



0(1), 
0(1), 



(3.52) 
(3.53) 



An+l) 
Ja,0 



A 



(n) 
0,0 



^■^0,0 ■ [|q,|_3,|q| + i],o 



0(ehl/2) + £(9(1) . o(er(|a|-3)/3l ) _ ^(^riai/si )^ 
0(el°l/2) + eO(l) • o(er(l"l-3)/3l ) _ o(en"l/3l )^ 



|a| ^ 3(n + 1), 

(3.54) 

\a\ > 3(n + 1), 



and 

An+l) 
JaA 



A 



(n) 



a.l 



^(n+1) I c-r("'+l) 

^■^0,1 ■ [|a|-3,|a| + l],l ^-^a.l " ^ [|a|-l,|a| + l],0 



o(ei«i/2+i) + ^o(£r(hhi)/3i ) + £0(^r(iai-i)/3i ) ^ o(^r(i"i+2)/3i )^ |„| ^3^ + 1^ 

0(el"l/2+i) + ^o(er(l-l-i)/3l ) + £o(^r(|ahi)/3l ) _ o(er(l"l+2)/3l )^ > 3n + 1. 

(3.55) 



This gives (lOHIl . The vahdation of (ITIil) and ([TIHD needs 



„(n+l) 
J a.k 



An) 
■1 a.k 



A 



[n) 
a,k 



A 



in) 

a,k 



A 



(n-l) 
a,k 



eiC 



(n+l) 
a.k 



V a,k 



A 



(n-l) 
a,k 



(n) 

[|a|-3,|a|+l],fc 
(n) 

[|a|-l,|Q!|+l],fc-l 
(n) 

[|o|-3,|a|+l],fc 



.(n+1) 



]o|-3,|a|+l],A;'' 



''a,k 



^a,k ' ^ [|a|-l,|Q|+l],fe-l/ 



[|a|-3,|Q|+l],fc'' 



siC 



(n+l) 



eL 



a,k 

(n) _ 

a,k 



,(n)N 
-a,k) 



.(n-l) 
[|a|-3,|a|+l],A:' 



(F 



(n) 

[|a|-l,|Q|+l],fe-l 



e(L 



(n+l) 
a,k 



-in). 
''a,k' 



F 



(n-l) 

|Q!|-l,|a|+l],A;-l' 



(3.56) 



and then (I?3U|1 . and P35D show that 

.(n+l) .(n) nz-^n+lN 



(3.57) 



The equations (IQSll — (fSlISD teh us that /^"^^ is never greater than 0(er(l"l+2fc)/3l ) for 
arbitrary n, and the leading order of fa,k appears at the [(|q;| + 2A;)/3]-th iteration step, 
and never changes later. Based on these results, we can remove some high order terms in 
the moment equations ()3.10p , and the remaining part is the formula for the regularization 
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term. Precisely, we have 



d=l 



^ fa-2ea,k + ^ RTtz 
d=l j=l 

+ k) + 2^ 2^ 77- (^^tr/, 



5t 



j=l d=l 



1 ^—-v d{RTt ) ^—-v 

+ 9 X] ^3^.' X] {RTtrfa-2ed-ej,k + Ujfa-2ea,k + (Oj + I) fa-2ea+e, ,k) 



(i=l 



k) ^ ga:"" (-^^tr/a-ej,fc-l + Ujfa,k-1 + (Oj + l)/a+ej,fc-l) 



-{Ga,k — fa,k) + /l-O.t. 



(3.58) 

where /i.o.t. stands for high order terms. The equations ()3.18p . ()3.19p and (|3.2ip can be 
used to make (|3.58p more compact. Noting that the terms containing fe^+e^fl or /g^.^i can 
be regarded as a high order term in ()3.8p , we can write p.58p as 



fa,k = £ 



E 



1 a(pRTa) 

p dxj 



fa— e,,fc RT^y 



dfa-ei,k\ 1 



3 3 



+ 5^^" ES^ E/« 



j=i / d=i 



^2ea,k 



^—^-g—^{RTtrfa-2ed-ej,k + (Oj + ^) fa-2e^+ej ,k) 
j=l d=l 



+ (m + /c) ^ ^'^^J™*'* {RTtrfa~e,,k~l + (Oj + l)fa+ej,k-l) 



i=l d=l 



If |a| = 2 and A; = 0, the above equation becomes 

3 



+ Gq,^ + h.o.t.. 



Suj duj 



+ (1 - Pr-^)/e.+e,,o + -6ijZ-'p{RT,^ - RTtr) + h.o.t.. 



With (13.81). it is reformulated as 



du/ 



(3.59) 



(3.60) 



(1 + <5i,)/e,+e,,o = -2Pr • epRTtr^^ + ()^^jPr • Z'^ p{RT,^ - RTtr) + /i.o.t., (3.61) 
which is similar as the Navier-Stokes law. Analogously, we can get the following relations 
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similar as the Fourier's law: 



Q^ = -tpRT,,^^§^ + h.o.t., 
fe„l = ^epRTtr + h.o.t. 

p.6ip and ()3.62p can be used to eliminate most partial derivatives in ()3.59p . Direct 
substitution gives the following result for \a\ ^ 3: 

3 



Ja,k — ^ / A a Ja-ei,k — -Kitr 7] 

^ 3 3 

+ . p X] + ^jd)fej+eafifa-ea-ej,k 

P j=l d=l 

1 ^ 
- -Z-\RT,^ - RTtr) fa-e,-e„k (3.63) 

i=i 

^ 3 3 

+ ^nT^T X] X] Qj[^'^tr:fa-2ea-e,,k + (Oj + I) fa-2ea+e, ,k] 
1 ^ 

+ urr Yl fe,,l[RTtrfa-e^,k-l + [Oj + l)/Q+e^,fc_l] + h.O.t.. 

Neglecting the high order terms and applying the result for |a| = Mq + 1, k = and 
\a\ = Ml + 1, k = 1, we carry out a reasonable approximation to the moments of the 
corresponding orders. This completes the closure of the moment system. 

There is still one question about the relation between Mq and Mi remaining. Since 
when \a\ ^ 1, fa,i always has the same order of magnitude as fa+ei+ej,Oi it is natural to 
choose Mo = Ml + 2. Thus the total number of moments is (Mq + 1)(M^ + 2Mo + 3)/3. 

3.4 Linearization of the regularization term 

The expression (j3.63p is rather complicated and is inconvenient for numerical imple- 
mentation. In [8j, we have used the technique of linearization to simplify the regularization 
term. Here, the same idea will be applied to derive a simplified regularization term. 

We consider a local problem where the distribution function is around the Maxwellian 
and the variations of the density, velocity and temperature are small. Thus we have the 
local expansions with a small parameter e indicating the magnitude of perturbation around 
the Maxwellian as 

p = po{l + ep), u = uo + e^/RTou, Tcq = ro(l + efcq), Ttr = To(l + eftr), , 

„ (3.64) 

x = Lex, p = Lpoy^oep, fa,k = Po{RTo)\''^/^+''eU,k for (a, fc) / (0, 0), 

where po, Uq and Tq are the reference density, velocity and temperature, respectively, 
the variables with hats " are dimensionless variables with magnitudes 0(1) and L is the 
characteristic length. The equations ()2.8p and ()3.20p show that 

p = poRTo{l + ep), Qj = poiRTof^hQj, j = 1,2,3, (3.65) 
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where p and Qj are also 0(1) dimensionless variables. Now we put ()3.64p and ()3.65p into 
(j3.63p . Eliminating the constants on both sides, we get 



- Pr(l + ep) ^ I^TT^ " + 



3 3 

e 



2Prfl + 6/51 ^^^^ ^ ^jd)fe,+ea,oL-ea- 
^ ' ' j=l d=l 



3 

fa—e^—Cj ,k 



3 3 

5(1 + e/))i?(l + eTtr) ^ ^ 

3 

(1 + ep)R{l + eTtr) ^ 



(3.66) 



Reserving only leading order terms on the right hand side, one has 



j=l J 

Using ()3.64p and p.65p again, we get a simple approximation of f^y- 

j=i J 

This approximation is to be applied to |a| = + 1 and used in our numerical method. 

Remark 1. The linearization may cause the loss of accuracy in the moment method. 
However, since the NRxx method is applicable up to arbitrary order of moments, the loss of 
accuracy can be got back by increasing the highest order in the system by 1. An important 
advantage of the regularization term (j3.68p is to smooth the profile of macroscopic variables 
such that the unphysical subshocks can be eliminated (see e.g. jl2l I20j ). We will find in 
the numerical experiments that the regularization term (|3.68p actually acts as a diffusion. 

3.5 Numerical method 

The framework of the numerical scheme for the NRxx method in the polyatomic case 
is generally the same as that in the monatomic case. The fractional step method is utilized 
to treat convection term and collision term separately. For the convection term, the finite 
volume method with the HLL numerical flux is employed. Below we consider the one- 
dimensional case and suppose a uniform spatial grid with grid size Ax is used. For an 
arbitrary quantity ■0, the symbol ij:'^ denotes the average value of •(/' on the i-th grid at the 
n-th time step. Then the whole algorithm is outlined as follows: 

1. Let n be zero and set /"(^,/) to be the initial value. 
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2. For each i, apply the technique of hnear reconstruction to determine the distributions 
on both boundaries of the i-th grid. The results are denoted as /j^/2(^'^) 

3. For each reconstructed distribution function, use ()3.68p for |a| = M^, + 1 to approx- 
imate the truncated moments. 

4. Use CFL condition to determine the time step At: 



At max ■ 



\u\^ + C{Mo)^{RTtr)f , 2(M + 1 



+ 



Ax ' (Ax)2 ^ I - 2 



{eRT,,)H ^-CFL, (3.69) 



where C{Mq) is the maximal root of the Hermite polynomial HeMo{x), and CFL is 
a specified Courant number between and 1. 

5. Solve the convection part with the HLL scheme: 



(3.70) 



where 



G'r+i/2(^) 



^ -^1+1/2' 



^r+1/2^1 /ill/2 (^'^) - ^r+i/2^i/m/2(^'-^) 



\n+ _ \n- 
^j+1/2 1+1/2 



+ ■ 



^r+l/2'^r+l/2[/rM-l/2(^'-^) /i+1/2 

it I)] 



1+1/2 t+1/2 



K+1/2 < < A"+^/2' 



^lfr+l/2itl), 



^ -^1+1/2' 



(3.71) 



and 



A"- = min {(ni):V~^/2 " C {M,) ^{RT^^^, {u,\ 



i+l/2 



n+ 
i+l/2 



X 



= max{(ni)^-,/2 + C{M^)^[{BT~fl^^, k^{flti/2 + ^^(^o)^^^?^^!^} • 

(3.72) 

6. Solve the collision-only equation by At, and the result is denoted as fl^~^^{^,I). 

7. Increase n by 1 and return to Step [2j 

The details of linear reconstruction and the numerical approximation of the truncated 
moments are almost the same as the monatomic case, and we refer the readers to [8] 
for our implementation. Here only Step [S] and Step [U] are expanded in the following 
subsections. 



3.5.1 The convection step 

In [7], we have mentioned that two operations are needed to accomplish the HLL 
scheme. One is the calculation of £,jf for / represented by (j3.5p . and the other is the 
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linear operation on distribution functions. The former is straightforward: 

aGN3 km 

= 2Z XI •^"'^ [RTtr'll^a+ej,k,Ttr,Ti^ti'", J) + Uj'lpa,k,Tu,Tir.ti'" ^ J) + aj1pa-e,,k,Ttr,T,^ti'"^ J)] ' 
oeN3 km 

(3.73) 

where the recursion relation of Hermite polynomials is used, and 

v = {i- u)/^/mr, J = I^'^/{RT,^t). (3.74) 

In order to make linear operations on distribution functions applicable, we have to 
solve the following problem: 

For f{^,I) represented by (|3.5p . find a series of coefficients f'^ ^, such that 



int , 



QGN3 km 

for some given u' , T/^, and Tj'^^. 
This can be tackled by two steps. First, we will find a series of coefficients Z^'^, such that 



fi^J) - Yi fa,k'^»,k,Ttr,T(^^ { ^/pr^' UT' I 

am-^km \VR^tv ^^intj 



(3.75) 



Note that only fa,kS with A; = 0, 1 are interested, which significantly reduces the difficulty 
of the problem. We calculate 



'a,k 



JR3xR+ 



I)^a,k,Ttr 



exp 



■ + 



1RT%r RTy 



ref 



d^dl 



(3.76) 

using both ()3.5p and ()3.75p . where T^ef is an arbitrary positive real number. Thanks to 
the orthogonality of Hermite and Laguerre polynomials, the results can be obtained as 



^a,0 — Cofa,0 — Cof'^Q, 



where 



/q!,1 + lj^{RTrci — RTmt)fa,0 



Ci 



fa,l + 2^R'^rei — RT(^t)fa,0 



Co = T7:747T7r(m+ l)(i?rtr)-(l"l+=^)(i?rref)-''/2 



5 (27r)3/2 

Ci = |_^r(m + 2)(i?Tt,)-(l"l+3)(i?T,ef)-(^/2+2)^ 



6 (2^)3/2 

With (f3T7l) and (fSlTS]) . we immediately get 



fa,0 — fa,0, fa,l — fa,l + ^(^^nt ~ RTmt)fa,0- 



(3.77) 

(3.78) 

(3.79) 
(3.80) 

(3.81) 
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The second step is to calculate ^ from /^'^. Since the scaling factor on the /-axis RT(^^^ 
is not changed in such transformation, we can use the very technique introduced in to 
obtain ^ efficiently. Specifically speaking, we define 

F{v, J,t)=Y,Y1 P-,kir)[iT - 1)t + l]l"l+=^V'a,fc,r...,7;;, {[{f - l)r + 1]^; + rw, j) , 

aGN3 fcGN 

(3.82) 



where 



If we require 



/ Tfr U — u' , , , , 

V -'tr V^^tr 



— EE 0, Vr G [0, 1], and F„,fc(0) = Z^'^^, V(a, A;) G x N, (3.84) 

then it is easy to find -Fa,fc(l) = fc, V(a, /c) G x N. The next job is to write ^ in the 
form of ()3.5p . and then require each coefficient to be zero. Thus an ordinary differential 
system is obtained. The calculation of ^ is almost a repetition of that presented in [7], 
which is omitted here. The resulting ordinary differential equations are 



dr 



dFa,k 



(3.85) 



where 



Sir) = . (3.86) 

(r-i)T + i 

In our implementation, we solve (j3.85p using Runge-Kutta method until r = 1, and then 
Fa,k{^) = f'ak be obtained. The readers can find some properties of ()3.85p in [7]. 
Equations ()3.8ip and p.85p give a practical way to calculate ^ for fc = 0, 1. 

3.5.2 The collision step 

It remains to give a numerical method to solve the collision-only equation 

f = i(G-/) (3.87) 

by one time step. The corresponding moment equations can be extracted from ()3.10p by 
removing the terms arising in the convection term. This results 

~g^ + 2^ ''''' 2 — di — Z_^J'^-'^^d,k-{rn + k) — }a,k-l = -{(-^a,k - Ja,k)- 

d=l d=l 

(3.88) 

Let (a, k) = (0, 0) and (a, k) = (ej, 0) respectively, and one has 
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Thus the second term in ()3.88p vanishes. Let {a,k) = (2ej,0), j = 1,2,3 and sum the 
equations up, and we get 

^(^eq-Ttr) (3.90) 



dt 

Setting (a, k) = (0, 1) in ()3.88p . it becomes 

dt eZ 

Using ()2.7p . it is easy to obtain 



cq 



dt 



0. 



(3.91) 



(3.92) 



Equations ()3.89p and (|3.92p agree with the fact that the density, velocity and temperature 
are not changed by collision. Now (j3.88p can be rewritten as 



dfa,k 

dt 



G. 



a,k 



fa,k — Z (i?Tcq — RTi 



(l ' \ 

n fa-2e,,k - {m + k)fa,k-l 

V7^^ J 



(3.93) 

Setting (a, k) = (e^ + ej, 0) in ()3.17p and substituting the result into ()3.93p . one finds that 
the collision-only equation for /e.+ej,o 

dfei+Sjfi _ 1 



has also a simple form: 



fe,+ej,0, hj — 1,2,3, 



(3.94) 



dt ePr" 

which agrees with the settings in [6J. Now suppose we want to solve (j3.87p from to 
^n+i ^ For simplicity, for any quantity 'ip, we use V'", V'"'"''^ ^ind ip^^^/"^ to denote ip{t'^), 
'ip(t"'~^^) and ip (^(t" + t""*"^)), respectively. Then the following relations holds analytically 
for i,j = 1,2,3 and t G 

p(,t)^p^, uj{t)^u^, n^{t)^T:^, 



rint(t) = r,"q + (7i^,-T,;)exp 



eZ 

t-f^ 
eZ 



(3.95) 



ePr 

These are deduced from IK89\\ — (^3M\\ and (ISTMl) . Based on (K95\\ . G„,fc(t) can also 
be obtained since it can be observed from (j3.1ip and (j3.15p — (j3.17p that Ga,k are fully 
determined by the variables listed in (|3.95p . For other cases, meaning \a\ > 2 if A; = and 
|a| > if /c = 1, the Crank-Nicolson scheme is employed to give a numerical approximation 
of I^Ml- 



J a,k J a,k J- 

At ~ e 



G 



Jn.k ^ J I 



i-n+l 
n+l/2 _ J a,k 
a,k 



a,k 



z-^ fi2r"+i/2-i?r"+^/2 



3 fn+i I fn 
i J a~2ej,k Ja-2ej,k 

9 2^ 3. 



fu+l I fn 
J a,k-l J a,k~l 



{m + k) 



(3.96) 



where At = t""*"^ — t". Note that no linear system needs to be solved when applying ()3.96p . 



since when solving /^^^, the terms /^^2e- k fak-i have always been obtained, and 
then (|3.96p is simply a linear equation of f[ 



n+l 
a,k ■ 
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Remark 2. When Pr = Z = 1, one has T = RT^qld and Trd = T^q in ()2.2p . Obviously, 
in this case, the ES-BGK model reduces to the BGK model, and the above technique for 
processing the collision terms is still valid. 

4 Numerical examples 

In this section, two one-dimensional numerical examples are presented to validate our 
algorithm. The spatial variable x will be written in plain font as x. For both tests, the 
non-dimensional form of the Boltzmann equation (|2.1|) is used. Thus the gas constant R 
is chosen as 1, and the Knudsen number Kn, which is the ratio of the mean free path 
to the characteristic length, controls the rarefaction of the gas. Only the diatomic gas is 
considered below; therefore 5 equals to 2.0. The CFL number is set to be 0.95 in all runs. 

4.1 Shock tube test 

There have been a number of studies on using the moment method to solve shock tube 
problems. In [2^, the 13-moment case is carefully investigated and the numerical results 
in O [23] indicate that the theory for 13-moment case can also be applied to systems with 
almost any number of moments. And in [71[9], the shock tube problem in the monatomic 
case is calculated to show the convergence of the original NRxx method and its improved 
version when the number of moments increases. Here the similar settings are used the test 
the polyatomic NRxx method. The initial conditions are 



where pi = 7, pr = 1, and = = 1. Recalling that R = 1 and 6 = 2, we find the 
fluid states on both left and right sides are in equilibrium. As an artificial test, a simple 
expression for the viscosity coefficient p is chosen as 



Four additional parameters, including the Prandtl number Pr, the relaxation collision 
number Z, the Knudsen number Kn, and the maximal moment order Mq, need to be 
defined for this problem. Below, different combinations of these parameters are tested in 
our numerical examples to show different properties of the polyatomic NRxx method. In 
all the experiments, the computational domain is [—2, 2] and discretized using 400 uniform 
spatial grids. In the following subsections, all plots show the numerical results at t = 0.3. 

4.1.1 Convergence in the number of moments 

In this part, we set Pr = 0.72 and Z = 5, and test the behavior of solutions for several 
Knudsen numbers when Mq increases. In order to provide a reference solution, Mieussens' 
conservative discrete velocity model (CDVM) [16\ [10] is computed. For CDVM, we use 
the technique of dimension reduction to speed up the computation. That is, we define 




X < 
X > 



(4.1) 



p = Kn ■ Tc 



(4.2) 




(4.3) 
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and solve g, hi and h2 instead of the distributions with full three-dimensional microscopic 
velocity. The velocity space is discretized by 400 grids. Currently, this skill is not used in 
the NRxx method. 

First, a relatively dense case Kn = 0.05 is considered. The density and equilibrium 
temperature results for Mq = 3, • • • ,8 are given in Figure [TJ For Mq ^ 5, the NRxx results 
match with the CDVM results very well. This agrees with the observation in [71 [9] that 
for small Knudsen numbers, a small number of moments can describe the macroscopic 
quantities in a high accuracy. The results for a rarefied case Kn = 0.5 from Mq = 3 
to Mq = 20 are presented in Figure [21 It is obvious that in order to match the CDVM 
results, much greater number of moments are needed. However, we can still find the NRxx 
results converge as Mq increases, and the limit is probably the solution of the Boltzmann 
equation. 

Now we consider a severe case Kn = 5, and the results are given in Figure [3l Although 
the NRxx solutions deviate from the reference solution greater than those in Figure [21 the 
convergence is again very clear. From Figure [T]—[3l we can find the theory in [23 is also 
valid for the regularized moment methods. For a fixed choice of Mq, the corresponding 
moment system always fails to describe the physical phenomenon when t — )• (or Kn — )• oo 
for a fixed time t) due to the very strong non-equilibrium. As t increases, the collision term 
starts to show an effect of dissipation, and the solution of the moment system gradually 
presents its physical meaning. As is shown in [3 |23], for a greater Mq, such progress is 
faster. Our numerical results correctly exhibit this behavior. 

4.1.2 Comparison between BGK and ES-BGK collision terms 

As is known, for monatomic gases, the BGK model fails to predict the correct Prandtl 
number, while Pr is considered as a parameter in the ES-BGK collision term. For poly- 
atomic gases, besides the Prandtl number, the BGK model also gives incorrect relaxation 
collision number Z. Actually, the BGK model always gives Z = 1, which means the 
translational and internal temperatures tend to the equilibrium temperature more rapidly 
than the ES-BGK model. Thus it can be expected that the BGK model gives incorrect 
translational temperature, internal temperature, and heat fiuxes. 

As a test, we set Pr = 0.72, Z = 5, Kn = 0.05 and Mq = 5, and both BGK and 
ES-BGK collision models are computed. The results are shown in Figure [H and Figure El 
In Figure [H it is found that the BGK model gives a pretty good prediction of the density, 
whereas the deviation of temperature between two models is significant. Figure [SI shows 
that the BGK model provides much smaller difference between the translational temper- 
ature and the internal temperature than the ES-BGK model, which indicates different 
relaxation collision numbers involved in the two models. The heat flux qi is defined as 



The difference in the heat flux is caused by the discordance of both the Prandtl number 
and the relaxation collision number. 

4.1.3 Comparison between the monatomic case and the polyatomic case 

Let Z = CO and define 




(4.4) 




(4.5) 
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X X 

(e) Mo = 7, 176 moments (f) Mo = 8, 249 moments 

Figure 1: Numerical results for the shock tube problem with Kn = 0.05. The dashed lines 
are the NRxx results, and the solid thin lines are the CDVM results with linearization. 
The dashdot lines are the results of discrete velocity model. The black lines denote the 
density p and the gray lines denote the equilibrium temperature Teq. 

Integrating the both sides of (|2.1|) over with respect to /, it is not difficult to find 
that the reduced distribution function g satisfies the monatomic Boltzmann equation with 
ES-BGK collision operator. Thus, it is natural to expect that when the relaxation collision 
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X X 

(e) Mo = 7, 176 moments (f) Mo = 8, 249 moments 



Figure 2: Numerical results for the shock tube problem with Kn = 0.5. The dashed lines 
are the NRxx results, and the solid thin lines are the CDVM results with linearization. 
The dashdot lines are the results of discrete velocity model. The black lines denote the 
density p and the gray lines denote the equilibrium temperature Tgq (continued on the 
next page). 

number Z gets greater, the polyatomic case will get closer to the monatomic case. The 
part is devoted to the numerical validation of this behavior. 
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X X 

(k) Mo = 13, 924 moments (1) Mo = 14, 1135 moments 

Figure 2 (continued): Numerical results for the shock tube problem with Kn = 0.5. The 
dashed lines are the NRxx results, and the solid thin lines are the CDVM results with 
linearization. The dashdot lines are the results of discrete velocity model. The black lines 
denote the density p and the gray lines denote the equilibrium temperature Teq (continued 
on the next page). 

The NRxx method for monatomic gases has been introduced in [3 [9l [8], where a 
BGK collision model is considered. For the monatomic ES-BGK model, the collision 
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X X 

(q) Mo = 19, 2680 moments (r) Mo = 20, 3101 moments 

Figure 2 (continued): Numerical results for the shock tube problem with Kn = 0.5. The 
dashed lines are the NRxx results, and the solid thin lines are the CDVM results with 
linearization. The dashdot lines are the results of discrete velocity model. The black lines 
denote the density p and the gray lines denote the equilibrium temperature Teq. 

only equation can be analytically solved, and the result will be reported elsewhere. Four 
relaxation collision numbers Z = 1, 10, 100, 1000 are considered here, and other parameters 
are Pr = 2/3, Kn = 0.01, Mq = 5. The numerical results can be found in Figure [6l It 



24 




-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 -1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 



X X 

(c) Mo = 9, 340 moments (d) Mo = 12, 741 moments 

Figure 3: Numerical results for the shock tube problem with Kn = 5.0. The dashed lines 
are the NRxx results, and the solid thin lines are the CDVM results with linearization. 
The dashdot lines are the results of discrete velocity model. The black lines denote the 
density p and the gray lines denote the equilibrium temperature Teq (continued on the 
next page). 



clearly shows that the polyatomic result tends to the monatomic result gradually as Z 
increases. 

4.2 Shock structure of nitrogen 

In this section, we will use the polyatomic NRxx method to compute the shock struc- 
ture of nitrogen, trying to reproduce the experimental results reported in [2]. In order to 
get a steady shock structure with Mach number Ma, we solve a Riemann problem with 
the following initial condition until a steady state: 



f m -r /t n = / Pl^^^^'TuTi ((^ - ui)/^/Ti, I/Ti) , x<0, 
J{u,x,!^,i) \ p^^oo,T,.,T,. ((^-^^r)/^/T;,//T.), :e > 0, 



(4.6) 
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X X 

(g) Mo = 21, 3564 moments (h) Mo = 24, 5225 moments 

Figure 3 (continued): Numerical results for the shock tube problem with Kn = 5.0. The 
dashed lines are the NRxx results, and the solid thin lines are the CDVM results with 
linearization. The dashdot lines are the results of discrete velocity model. The black lines 
denote the density p and the gray lines denote the equilibrium temperature Teq. 

where 

pi = l, = (V7^«,0,0f , Ti = l, 

_ (7 + l)Ma^ 2jMa^ - (7 - 1) (^•'^) 

(7-l)Ma2 + 2' "'■-p/'' (7 + l)p. ■ 

Here 7 is the adiabatic index. For nitrogen, 7 equals to 1.4. The Prandtl number is chosen 
as 0.72 as in [24J. The Knudsen number is Kn = 0.1, and the grid size is Ax = 0.005. 
The computational domain is [—1.5, 1.5], which is large enough to cover the whole shock 
structure. 

It remains to give the expressions of Z and fi. They have significant influence on the 
thickness of the shock. For the relaxation collision number Z, both the gas-kinetic model 
j24j and the direct simulation Monte Carlo (DSMC) [2] show that Z = 4 01 Z = 5 best 
fits the experimental data. However, in both [13j and [3j, where the Rykov and ES-BGK 
models are used respectively, it is reported that a smaller Z between 2 and 3 gives better 
numerical results. The same conclusion is drawn by our numerical experiments. Since [13] 
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Ql 1 < < < < < 1 < 1 ^0.7 

-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 

X 

Figure 4: Comparison between BGK and ES-BGK models. The black lines are the results 
of the ES-BGK model, and the gray lines are the results of the BGK model. The solid 
lines denote the profile of density p, and the dashed lines denote the profile of equilibrium 
temperature T^q. 




-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 

X 

Figure 5: Comparison between BGK and ES-BGK models. The black lines are the re- 
sults of the ES-BGK model, and the gray lines are the results of the BGK model. The 
dashed and dotted lines denote the profile of translational temperature Ttr and internal 
temperature Tint, respectively, and the solid lines denote the heat flux qi. 



also considers the shock structure problem, we use the same settings here: 
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-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 



(a) Density profile 




(b) Translational temperature profile 
Figure 6: Comparison between monatomic and polyatomic cases. 



Six Mach numbers ranging from Ma = 1.53 to Ma = 6.1 are taken into account. 
Similar as [S], in order to avoid the problem of hyperbolicity, only the 24 moment system 
(Mq = 3) is used in our computation. The numerical results are plotted in Figure [71 where 
all macroscopic variables are normalized so that the computational results can match the 
data in [2J. Precisely, we use 



P 



P- Pi 



Ti 



'-int 



Ti 



and A denotes the mean free path. It can be found that the density profiles are in very 



(4.9) 
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good agreements with the experimental data and no subshocks exist in the shock struc- 
ture. With increasing Mach number, the numerical result gradually deviates from the 
experimental data. Only when Ma is as great as 6.1, the deviation in the low density 
region (around x/A G (—4, —1) in the figure) is becoming obvious. 



o p. exp. data 



- - ■ T„ 

o p. exp. data 
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5 Concluding remarks 



In this paper, the NRxx method is extended to the polyatomic gases. Further inves- 
tigations such as the boundary conditions and the multidimensional simulations are in 
progress. 
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A Properties of Hermite and Laguerre polynomials 

The Hermite polynomials defined in ()3.4p are a set of orthogonal polynomials over the 
domain (—00, -|-cxd). Below we list some of their properties which are used in this paper: 

1. Orthogonality: Jj^ i/e„,^(x)i/e„2(x) exp(— x^/2) dx = ni!\/27r(5nin2; 

2. Recursion relation: i?e„+i(x) = xife„(x) — nife„_i(x); 

3. Differential relation: He'^{x) = nHen~i{x). 

All these properties can be found in many mathematical handbooks such as [T]. And the 
following equality can be derived from the last two relations: 

[Hen{x) exp(-xV2)]' = -Hen+i{x) exp(-xV2). (A.l) 

As introduced in Section [3. H the Laguerre polynomials defined in (|3.3|) are orthogonal 
over [0,-|-oo). The Laguerre polynomials are closely related to the Hermite polynomials, 
and they have very similar properties: 

1. Orthogonality: l[.™^(x)L^™^(x)x"' exp(-x) dx = T^'^^^ifca; 

2. Recursion relation: {k + 1)l[,™^-^(x) = (m + 1 + k — x)l[,'"^(x) — xL^j^^^\x); 

3. Differential relation: [L^^\x)]' = —L^^'^^\x). 
And the last two relations give 

[4'"^(x) exp(-x)]' = x-^iik + 1)471(2^) -im + l + A:)4'"^(^)] exp(-x). (A.2) 

B The deduction of polyatomic moment equations 

In this appendix, we are going to give the detailed deduction of (|3.10p . For simplicity, 
we define 

^3 1+3 ^ / v'^\ 

^l,a,TM = (V2^) (RTtr)-^ n H^c^M exp ( ) , 

d=i ^ ^ (B.l) 
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Thus 'ipa,k,Ttr,Tinti'"^ = 'ipi,a,Ttr{v)'ip2,k,Tint{J)- It has been deduced in [8] that 



E 



dud I 



where 77 stands for t or Xj, j = 1, 2, 3. Now using ()A.2p . we have 



(B.2) 



drj 



(m) 



RTmt j 



j2/5 



{k + i)Li':i 



Since m = 5/2 — 1 and 



j2/5 



(m + 1 + fc)Ll'") 



exp ^- 

/ j2/5 



\ RTinX, 
RT\nt I 



exp 



RT; 



int 



int 



int 



P/& drj \RTi^t^ 
the equation ()B.3p can be simphfied as 



1 g(i?rmt) 



^ / 



(m) 

7fc 



j2/S 
RTint 

-1 



j2/5 



RTi 



exp 



int 



int 



Noting that 



7^+} _ r(m + k + 2) T{k + 1) _m + k + l 



r(m + /c + l) r(A; + 2) k + l 
one finally obtains a simple expression: 



^ / 

9;^V'2,.,T.„ 



j2/S 



[m + k -\- 1)^/'; 



j2/5 



^^int / 

Thus the derivative of the basis function (j3.ip is 
d 



'2,fc + l,Tin 



int 



(i=l 



-^^Q+ed,A;,Tt„Tint + 2 'q:^ Va+2ea,k,Ttr,Ti„ 



(B.3) 



(B.4) 



(B.5) 



(B.6) 



(B.7) 
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Here ipa,k,Ttr:,Tint stands for 



(B.9) 



int 



The parameters are omitted for conciseness. 

Now we expand the left hand side of (|2.ip into series. Using (jB.Sp . one immediately 

has 



dt 

= EE 

agN3 km 



dfa,k , , f ^ I 

—g^Va,k,Ttr,Tint + Ja,k-Q^Va,k,Ttr,Ti„ 



dfa,k sr^ dUd , 

dt dt ^'"^'"'^ 

d=l 

+ 2^ Ja-2ed,k - [m + k) — 



(B.IO) 



2 dt 



dt 



1pa,k,Ttr,Ti^ 



For the convection term, we have 



j=l ^ j=l aeNS km 

3 



dfa,k sr^ Oud 

d=l ■' 



+ ^ 2^ Ja-2ea,k - [m + k) — 



2 dx 



^ d=l 



^a,k,Ttr,Ti^t- 

(B.ll) 



Now we use ()3.73p and get 

3 



aGN3 km j=l 

dud 



Tjrr '^fa-ej,k dfa,k . , . ^^'^fa+e.j,k 



d=l ^ 

+ {RTtrfa-2ea-ej,k + Ujfa-2ea,k + ("j + l)/Q_2ed+ej 



2 dx 
(m + /c) 



d=l 

dxj 



[RTtvfa~ej,k-l + Ujfa,k-1 + (Oj + Ij/o+ej.fc-l 



(B.12) 



Cohecting (jRlOl) (jRT2|) and (fS^Q]) . the moment system (f3J0]l fohows naturally. 



C Expansion of the generalized Gaussian 



This section is devoted to the calculation of Ga,o, which is defined in (|3.9p . In this 
appendix, G is considered as a function of $, and /, where the parameters t and x are 
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omitted. It can be deduced from the orthogonality of Hermite and Laguerre polynomials 
that 

Ga,0 = Ca,T,. [ Pa,TMG (u + y^rV, l) dv dl , (C.l) 

where Pa,Ttr is a polynomial defined as 

) = (V2^) (RT,,)-^ H He^M, (C.2) 
^ d=l 



and Ca,Ttr is a constant dependent on a and Tt^: 



(2vr)-|(CTt,)3+l^l 

L^a.Ttr = j j i • l^-^j 

For i G {1, 2, 3}, if Oj > 0, the recursion relation of Hermite polynomials shows that 

Pa,TM = iRTtr)--2ViPa-e,,TM - {RT,,)-\ai - lK_2e„T,.(^). (C.4) 

Noting that 

ai[ak-Oik) ai[ai-l) 
one directly obtains from (|C.ip that 

G„,0 = C„,T,,(i?Tto)-^ / ViPa-e^^rMG (u + ^/RT[,V,i) dvdl - ^G^-2e„0. 

(C.6) 

Now the expression of G (|2.2p is put into the above equation. After integrating with 
respect to /, one has 

'^"'0 = / , ^rn / 'JiP»-e.,nAV) exp —V T V] dv Ga-2e,,0- 

^ydet{2^^T) Jr^ \ I J at 

(C.7) 

The matrix T is required to be positive definite, since the density of the fluid should be 
finite. Thus, there exists a matrix TZ = {vij) such that 7~ = {RTt^)Tni'^. Making the 
transformation w = W^v, and noting that det(r) = (i^Ttr)^[det(7^)]^ ^f]7!\ becomes 

Ga,Q = TtS^TF^^^^'^ I ^jP<^-^^,T,AT^w) exp(-|'U;|V2) dw - ^^Ga-2e,fi. 

[V27ry{RTtry jri Jr^ ai 

(C.8) 

The following relation is a direct result of the differential relation of Hermite polynomials: 
d 



dw 



-pa,T,A'^w) = (RTtr)-^ Yl akTkjPa-e.^T.AT^w). (C.9) 
^ k=l 



Thus it can be obtained by integrating by parts that 

/ WjPa-e„Ttr{T^w)exp{-\w\'^/2)dw 
Jrs 

3 (C.IO) 

= (i?rtr)~5 V(afc - 6ik)rkj / Pa-ei-ek^TtX'^w) exp(-|u;| V2) dw. 



k=l 
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Now we substitute (|C.10p into (jC.Sp . and apply the transformation v = TZw. The result 
is 

■ v^^^)U h ' (c.ii) 



/ Pa-e,-efc,Tt.(^^)exp —V T V \ dv G, 

Using (jC.Sp . it is not difficult to find 



a-2ei,0- 



Gafl = j ^ ^jj ^ ^feiGcf-e.-efe — Ga-2ei,0 j • (C-12) 

\j=l k=l 



Recalling T = {RTx.r)TZTZ , the above equation can be written as 

3 

Oii 



1 ^ 

Gafl = ^(Ajfc — RTtrSik)Ga-ei-ekfl, (C.13) 



k=l 



where Xik is the (i, A;)-element of T- The final result p.l7p is then obtained by substituting 
the detailed expression of T into ()C.13p . 
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